This script merges CalMAPPER activity data with treatment polygons. This is a necessary step before analyzing prescribed fire data in CalMAPPER.
calmapper_dir <- fs::dir_ls(ref_path, recurse = T, glob = '*CalMAPPER', type = 'directory')
calmapper_gdb <- fs::dir_ls(calmapper_dir, recurse = T, glob = '*.gdb')
if(length(calmapper_gdb) > 1){
stop( "There's more than one CalMAPPER .gdb file. script assumes only 1.")
}st_layers(calmapper_gdb)## Driver: OpenFileGDB
## Available layers:
## layer_name geometry_type features fields
## 1 CMDash_ProjectTreatments Multi Polygon 1560 21
## 2 CMDash_TreatmentPols Multi Polygon 2784 23
## 3 CMDash_TreatmentLines Multi Line String 23 22
## 4 CMDash_TreatmentPnts Multi Point 136 20
## 5 CMDash_Activities NA 20439 38
## 6 CMDash_Metadata NA 1 7
## 7 CalMapper_fire_act Multi Polygon 2799 28
## crs_name
## 1 WGS 84 / Pseudo-Mercator
## 2 WGS 84 / Pseudo-Mercator
## 3 WGS 84 / Pseudo-Mercator
## 4 WGS 84 / Pseudo-Mercator
## 5 <NA>
## 6 <NA>
## 7 WGS 84 / Pseudo-Mercatoract <- st_read(calmapper_gdb, layer = 'CMDash_Activities')## Reading layer `CMDash_Activities' from data source
## `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CALFIRE spatial data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb'
## using driver `OpenFileGDB'## Warning: no simple feature geometries present: returning a data.frame or tbl_dftrt <- st_read(calmapper_gdb, layer = 'CMDash_TreatmentPols')## Reading layer `CMDash_TreatmentPols' from data source
## `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CALFIRE spatial data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb'
## using driver `OpenFileGDB'## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts. The
## processing may be really slow. You can skip the processing by setting
## METHOD=SKIP, or only make it analyze counter-clock wise parts by setting
## METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock
## wise defined## Simple feature collection with 2784 features and 23 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13842530 ymin: 3842624 xmax: -12941530 ymax: 5161286
## Projected CRS: WGS 84 / Pseudo-Mercatoract <- act %>%
filter(!ACTIVITY_DESCRIPTION %in% c("GIS Validation", "Education Outreach", "Project Administration", "Planning Meeting", "Public Contacts", "Water Site Development"))act %>%
group_by(ACTIVITY_DESCRIPTION, UNIT_OF_MEASURE) %>%
count() %>%
arrange(desc(UNIT_OF_MEASURE)) %>%
print(n=47)## # A tibble: 47 × 3
## # Groups: ACTIVITY_DESCRIPTION, UNIT_OF_MEASURE [47]
## ACTIVITY_DESCRIPTION UNIT_OF_MEASURE n
## <chr> <chr> <int>
## 1 Biomass Removal (Bone Dry Tons) Tons 85
## 2 Dozer Line Miles 35
## 3 Hand Line Miles 108
## 4 Road Grading Miles 6
## 5 Air Curtain Burner Hours 13
## 6 Boundary Mapping Hours 43
## 7 Environmental Review Hours 110
## 8 Limbing and Bucking Hours 1204
## 9 Public Meetings Hours 21
## 10 RPF Supervision Hours 86
## 11 Site Assessment Hours 93
## 12 Site Preparation (Manual) Hours 77
## 13 Site Preparation (Mechanical) Hours 7
## 14 Site Preparation (RxBurn) Hours 187
## 15 Trees Felled (> 6in dbh) Each 1026
## 16 Broadcast Burn Acres 868
## 17 Chaining Acres 100
## 18 Chipping Acres 2216
## 19 Commercial Thinning (Cable Yarding) Acres 20
## 20 Commercial Thinning (Tractor Yarding) Acres 42
## 21 Crushing Acres 51
## 22 Cultural Burning Acres 5
## 23 Erosion Control Acres 12
## 24 Follow up - Herbicide Acres 58
## 25 Follow up - Other Acres 6
## 26 Follow up - Slash disposal Acres 177
## 27 Fuel Break (Shaded) Acres 55
## 28 Grazing Acres 52
## 29 Herbicide (Post-Treatment) Acres 43
## 30 Herbicide (Pre-Treatment) Acres 17
## 31 Invasive Plant Removal Acres 7
## 32 Land Conservation Acres 4
## 33 Lop and Scatter Acres 717
## 34 Mastication Acres 1243
## 35 Pile Burning Acres 1926
## 36 Piling (Manual) Acres 2533
## 37 Piling (Mechanical) Acres 375
## 38 Pruning Acres 376
## 39 Rangeland Mowing Acres 61
## 40 Release - Herbicide Acres 16
## 41 Release - Mechanical Acres 99
## 42 Release - Other Acres 10
## 43 Site Preparation (CFIP) Acres 71
## 44 Thinning Acres 287
## 45 Thinning (Manual) Acres 4400
## 46 Thinning (Mechanical) Acres 442
## 47 <NA> <NA> 1act <- act %>%
select(-UNIT_OF_MEASURE, -PROJECT_STATUS, -QUANTITY, -ACTIVITY_STATUS, -GROUND_DISTURBING, -ThisFiscal, -LastFiscal, -PrevFiscals, -OBJECTID_1, -CFIP_FR, -CONTRACT)act_sf <- right_join(trt %>% select(PROJECT_ID, TREATMENT_ID),
act)## Joining with `by = join_by(PROJECT_ID, TREATMENT_ID)`act_sf %>%
st_drop_geometry() %>%
group_by(ACTIVITY_DESCRIPTION) %>%
count() %>%
print(n=50)## # A tibble: 47 × 2
## # Groups: ACTIVITY_DESCRIPTION [47]
## ACTIVITY_DESCRIPTION n
## <chr> <int>
## 1 Air Curtain Burner 13
## 2 Biomass Removal (Bone Dry Tons) 85
## 3 Boundary Mapping 43
## 4 Broadcast Burn 868
## 5 Chaining 100
## 6 Chipping 2216
## 7 Commercial Thinning (Cable Yarding) 20
## 8 Commercial Thinning (Tractor Yarding) 42
## 9 Crushing 51
## 10 Cultural Burning 5
## 11 Dozer Line 35
## 12 Environmental Review 110
## 13 Erosion Control 12
## 14 Follow up - Herbicide 58
## 15 Follow up - Other 6
## 16 Follow up - Slash disposal 177
## 17 Fuel Break (Shaded) 55
## 18 Grazing 52
## 19 Hand Line 108
## 20 Herbicide (Post-Treatment) 43
## 21 Herbicide (Pre-Treatment) 17
## 22 Invasive Plant Removal 7
## 23 Land Conservation 4
## 24 Limbing and Bucking 1204
## 25 Lop and Scatter 717
## 26 Mastication 1243
## 27 Pile Burning 1926
## 28 Piling (Manual) 2533
## 29 Piling (Mechanical) 375
## 30 Pruning 376
## 31 Public Meetings 21
## 32 RPF Supervision 86
## 33 Rangeland Mowing 61
## 34 Release - Herbicide 16
## 35 Release - Mechanical 99
## 36 Release - Other 10
## 37 Road Grading 6
## 38 Site Assessment 93
## 39 Site Preparation (CFIP) 71
## 40 Site Preparation (Manual) 77
## 41 Site Preparation (Mechanical) 7
## 42 Site Preparation (RxBurn) 187
## 43 Thinning 287
## 44 Thinning (Manual) 4400
## 45 Thinning (Mechanical) 442
## 46 Trees Felled (> 6in dbh) 1026
## 47 <NA> 1act_sf_fire <-
act_sf %>%
filter(ACTIVITY_DESCRIPTION %in% c("Broadcast Burn", "Cultural Burning", "Pile Burning"))test <- act_sf_fire %>%
mutate(ACTIVITY_START_test = as.Date(ACTIVITY_START)) %>%
mutate(ACTIVITY_END_test = as.Date(ACTIVITY_END)) This should be 2/11 - 2/11. The values that R automatically reads in from ArcGIS are one day off.
test %>%
filter(AQ_ID == 69461) %>%
select(ACTIVITY_START, ACTIVITY_START_test, ACTIVITY_END, ACTIVITY_END_test)## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13709070 ymin: 4748703 xmax: -13705400 ymax: 4751503
## Projected CRS: WGS 84 / Pseudo-Mercator
## ACTIVITY_START ACTIVITY_START_test ACTIVITY_END ACTIVITY_END_test
## 1 2020-02-10 16:00:00 2020-02-11 2020-02-10 16:00:00 2020-02-11
## SHAPE
## 1 MULTIPOLYGON (((-13706330 4...act_sf_fire <- act_sf_fire %>%
mutate(ACTIVITY_START = as.Date(ACTIVITY_START)) %>%
mutate(ACTIVITY_END = as.Date(ACTIVITY_END)) act_sf_fire <- act_sf_fire %>%
mutate(DURATION = ACTIVITY_END - ACTIVITY_START+1)act_sf_fire <- act_sf_fire %>%
mutate(YEAR = year(ACTIVITY_END))act_sf_fire %>%
filter(AQ_ID == 69461) %>%
select(ACTIVITY_START, ACTIVITY_END, DURATION, YEAR) ## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13709070 ymin: 4748703 xmax: -13705400 ymax: 4751503
## Projected CRS: WGS 84 / Pseudo-Mercator
## ACTIVITY_START ACTIVITY_END DURATION YEAR SHAPE
## 1 2020-02-11 2020-02-11 1 days 2020 MULTIPOLYGON (((-13706330 4...cm <- act_sf_fire
save(cm, file = "Rdata/CalMapper_activities_fire.Rdata")write.csv(act_sf_fire %>% st_drop_geometry(), file = "~/Reference data/CALFIRE spatial data/CalMapper/activities_fire.csv", row.names = F)
write.csv(trt %>% st_drop_geometry(), file = "~/Reference data/CALFIRE spatial data/CalMapper/treatments_fire.csv", row.names = F)act_broadcast <- act_sf_fire %>%
filter(ACTIVITY_DESCRIPTION == "Broadcast Burn")trt %>%
st_drop_geometry() %>%
group_by(TREATMENT_OBJECTIVE) %>%
count()## # A tibble: 5 × 2
## # Groups: TREATMENT_OBJECTIVE [5]
## TREATMENT_OBJECTIVE n
## <chr> <int>
## 1 Broadcast Burn 559
## 2 Forestland Stewardship 296
## 3 Fuel Break 215
## 4 Fuel Reduction 1542
## 5 Right of Way Clearance 172trt_broadcast_complete <- trt %>%
filter(TREATMENT_OBJECTIVE == "Broadcast Burn") %>%
filter(ACTIVITY_STATUS == "Complete")act_broadcast that
aren’t in trtcheck <- act_broadcast %>%
filter(!TREATMENT_ID %in% trt$TREATMENT_ID) %>%
nrow()check2 <- act_broadcast %>%
filter(!TREATMENT_NAME %in% trt$TREATMENT_NAME) %>%
nrow()if(check==0 & check2 == 0){
paste("There are no TREATMENT_ID values or TREATMENT_NAME values in activities that aren't in trt, at least for broadcast burns")
} else{
print("Alert! There are TREATMENT_ID or TREATMENT_NAME values in activities data that aren't in trt")
}## [1] "There are no TREATMENT_ID values or TREATMENT_NAME values in activities that aren't in trt, at least for broadcast burns"not_in_act <- trt_broadcast_complete %>%
filter(!TREATMENT_NAME %in% act_broadcast$TREATMENT_NAME) %>%
nrow()print(paste("Alert! There are", not_in_act, "records in trt_broadcast_complete that aren't in activities. These are potential missing data from the final output!"))## [1] "Alert! There are 7 records in trt_broadcast_complete that aren't in activities. These are potential missing data from the final output!"missing_from_activities <- trt_broadcast_complete %>%
filter(!TREATMENT_ID %in% act_broadcast$TREATMENT_ID) %>%
select(CALMAPPER_ID, PROJECT_NAME, TREATMENT_NAME, TREATMENT_OBJECTIVE, PROJECT_STATUS, ACTIVITY_STATUS, PROJECT_START_DATE, PROJECT_END_DATE, TREATMENTAREA_ACRES)
missing_from_activities %>%
st_drop_geometry() %>%
summarize(missing_acres = sum(TREATMENTAREA_ACRES))## missing_acres
## 1 653.99write_excel_csv(missing_from_activities %>% st_drop_geometry(), "calmapper_trt_not_act.xls")map <- mapview(list(trt_broadcast_complete, act_broadcast, missing_from_activities), col.regions=list("red","blue", "green"),col=list("red","blue", "green"))
map